On the effective permittivity of finite inhomogeneous objects 
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A generalization of the S-parameter retrieval method for finite three-dimensional inhomogeneous 
objects under arbitrary illumination and observation conditions is presented. The effective per- 
mittivity of such objects may be rigorously defined as a solution of a nonlinear inverse scattering 
problem. In this respect the problems of S-parameter retrieval, effective medium theory, and even 
the derivation of the macroscopic electrodynamics itself, turn out to be all mathematically equiv- 
alent. We confirm analytically and observe numerically effects that were previously reported in 
the one-dimensional strongly inhomogeneous slabs: the non-uniqueness of the effective permittivity 
and its dependence on the illumination and observation conditions, and the geometry of the object. 
Moreover, we show that, although the S-parameter retrieval of the effective permittivity is scale-free 
at the level of problem statement, the exact solution of this problem either does not exist or is not 
unique. Using the results from the spectral analysis we describe the set of values of the effective 
permittivity for which the scattering problem is ill-posed. Unfortunately, real nonpositive values, 
important for negative refraction and invisibility, belong to this set. We illustrate our conclusions 
using a numerical reduced-order inverse scattering algorithm specifically designed for the effective 
permittivity problem. 

PACS numbers: 41.20.Jb, 42.25.Bs, 78.20.Bh. 78.20.Ci 



I. INTRODUCTION 



The notion of an effective permittivity (conductivity, 
permeability) has been introduced within the effective- 
medium theory and was meant as a simplifying approx- 
imation of the scattering model for objects exhibiting 
inhomogeneity on scales much smaller than the scale of 
the spatial variation of the incident field Thus, a 
gas, a colloid, a suspension, or a powder mixture could 
be modelled as homogeneous media with some effective 
permittivities when interacting with a field having wave- 
length much larger than the size of the constituent par- 
ticles and the average distance between them. Extension 
of the effective medium approximation beyond its nat- 
ural limits of applicability, i.e. for smaller wavelengths 
where the multiple interparticle scattering is more pro- 
nounced, leads to effective parameters exhibiting strange 
and sometimes exotic behavior. For example, if we insist 
on describing the Earth atmosphere as a homogeneous 
medium, then the naive explanation of the sky's blue 
color could be a particular frequency-dependence of this 
homogeneous effective permittivity, so that it filters out 
other frequencies of the visual spectrum. In reality, as we 
know, the opposite happens. The high-frequency blue 
light is Rayleigh-scattered by the particles (molecules) 
of the atmosphere, so that even the parts away from 
the direct optical path to the Sun shine with blue light 
that reaches our eyes. Whereas the low-frequency red 
light passes the atmosphere almost without scattering 
and causes the red color of the sunset. This illustrates an- 
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other important property of all effective scattering mod- 
els. Should we decide to view the sky as a homogeneous 
medium, then we would have to introduce at least two 
effective models - one for the day and another for the 
evening, i.e. our effective model would depend on the illu- 
mination and observation conditions. It was also noticed 
in [2j that the effective model of an electrically inhomoge- 
neous atmosphere may show non-zero magnetic contrast 
- a conclusion of great importance to the present studies 
of composite media. 

In recent years the effective medium approach has been 
extensively used beyond its natural applicability limits, 
especially for modelling the response of metamaterials 
and photonic crystals. Moreover, the effective param- 
eters of such composite media seem nowadays to carry 
more physical meaning than a mere simplifying approx- 
imation of the scattering model. For example, the fol- 
lowing practical question may be asked: can a composite 
medium be used in the same way as one would use, say, 
a borosilicate glass, to cut out a lens or a light-bending 
coating (invisibility cloak)? Or is there something about 
strongly inhomogeneous composites that makes them dif- 
ferent from ordinary dense media in this respect? It 
should be mentioned, of course, that the very notions of 
the macroscopic permittivity and permeability of those 
"ordinary" media are themselves the product of the so- 
called continuum approximation, which is essentially an 
effective medium theory applied within its "natural" lim- 
its. Hence, the problem is, in fact, older and more funda- 
mental. It concerns the applicability, accuracy, and the 
physical meaning of the effective modelling as such. 

Partly due to the popularity of metamaterials, the 
derivation techniques and even the very concept of effec- 
tive permittivity (permeability) have been recently revis- 
ited by many authors [3-32]. There are two main ways 
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in which one can derive the effective permittivity of an 
inhomogeneous object. One is the traditional approach 
of macroscopic electrodynamics known as homogeniza- 
tion. It attempts to derive a local constitutive relation 
between the averaged field and a simplified permittiv- 
ity (conductivity, permeability) function, e.g. averaged 
over a representative cell or simply constant. Applica- 
tion of the homogenization to metamaterials is reviewed 
in [lol [3l| . The authors recognize the problem, which 
the relatively large and strongly scattering metamaterial 
particles pose for the applicability of this traditional ef- 
fective medium approach and view homogenization as a 
complementary method to the so-called S-parameter re- 
trieval technique. The latter technique, whose applica- 
tion to metamaterials was pioneered in ||, is the other 
popular way of looking at the problem where instead of 
explicitly deriving a local constitutive relation one is sim- 
ply matching the observed (simulated or measured) field 
from a composite inhomogeneous slab to the field from 
a homogeneous slab of the same thickness. Often this is 
also the way the exotic properties of metamaterials and 
photonic slabs are verified experimentally. From the the- 
oretical point of view the S-parameter retrieval method 
is very attractive, because it seems to be immune against 
the aforementioned scalability problems inherent to ho- 
mogenization techniques. Indeed, one does not care what 
the relative scale of inhomogeneities is, as long as there 
is a match with a field from a homogeneous object of the 
same outer shape. There exists a third method, which 
includes the measured internal field in the matching pro- 
cedure However, as we show here, it is essentially an 
extension of the S-parameter retrieval technique. 

Despite its generality, the S-parameter retrieval poses 
other questions. First of all, the retrieved effective pa- 
rameters for a homogeneous slab turn out to be non- 
unique [1, 0] . Secondly, their values are sensitive to the 
location of the slab boundary Q , orientation and regular- 
ity of the cells Q, and depend on the angle of incidence 
of the illuminating plane wave Q . This dependence on 
the wavevector k is considered by some to be the sign 
of anisotropy and spatial dispersion [2l| and has lead 
others [H, |23| . [3l[ to question the usefulness of the ef- 
fective medium parameters as such, since they do not 
serve their original purpose of simplifying the propaga- 
tion model any more. Also, the complex effective permit- 
tivity may sometimes show a negative loss, i.e. gain, [J] 
and larger than expected (from homogenization) positive 
losses (3. fllj] . 

The theoretical analysis of the S-parameter retrieval 
method has been limited so far to the slab-like configu- 
rations where an analytical solution for a homogeneous 
case is readily available. Here we consider a generaliza- 
tion of this approach to finite three-dimensional struc- 
tures, inhomogeneous, and not necessarily metamateri- 
als. We apply the volume integral formulation and show 
that we deal with a special kind of inverse scattering 
problem. We demonstrate the mathematical equivalence 
of the generalized S-parameter retrieval (or effective in- 



version) problem to the original problem of the effective 
medium theory. Despite the lack of explicit analytical 
solutions in 3D, a combination of spectral analysis and 
inverse scattering theory confirms all the anomalous fea- 
tures of the S-parameter retrieval method in the present 
general 3D case and sheds new light on their mathemat- 
ical origins. 

The most important and somewhat surprising results 
of our study are: in general, an exact effective model of 
lower complexity does not exist; if an effective permittiv- 
ity exists, it is not unique, and there is often no way one 
can choose a particular value on "physical" grounds; the 
goals of having an accurate effective scattering model and 
a unique effective permittivity contradict each other; the 
effective model is singular for (infinitely) many real non- 
negative values of permittivity (at least for electrically 
small objects). In most of the paper we limit ourselves 
to the retrieval of an effective permittivity of a homoge- 
neous (lower complexity) effective model, although, effec- 
tive models of higher complexity (anisotropic, magnetic) 
are briefly discussed as well. We illustrate the key effects 
on a number of numerical examples using an algorithm 
specifically designed for the problem. This algorithm is 
of interest in itself as it eases the computational bur- 
den of having to solve many forward scattering problems 
(one for each trial value of the effective permittivity) by 
employing the shift invariance of the Arnoldi iterative 
scheme. 



II. ANALYSIS OF THE FORWARD 
SCATTERING PROBLEM 

A. The volume integral equation method 

The most general definition of an effective scatterer 
is, in fact, quite simple and can be formulated as follows. 
An effective scatterer has the same or similar outer shape 
as the original object; different, possibly simpler, inter- 
nal structure; and produces the same (or approximately 
the same) scattered field as the original under identical 
illumination and observation conditions. In the following 
section we shall make this definition more precise. For 
the moment, let us consider a general inhomogeneous, 
isotropic, finite object occupying the spatial domain D. 
Its constitutive parameters are the spatially varying com- 
plex, possibly dispersive, dielectric permittivity e(x, w) 
and constant magnetic permeability fig. The surround- 
ing medium is infinite and homogeneous with the vacuum 
parameters £o and /io- Thus, the object has no magnetic 
contrast with respect to the background, whereas its elec- 
tric contrast is given by the function 



Obviously, x(x, u>) = 0, if x ^ D. Let the object be 
illuminated by an external source, which generates the 
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known electric field E m (x, ui) in vacuum. The total elec- 
tric field E(x, uj) inside the scatterer satisfies the follow- 
ing volume integral equation: 

E(x,w)- / G(x-x',w)x(x',w)E(x',w)dx' 

An (2) 
= E in (x,w), xefl, 

where G(x, uj) denotes the Green's tensor. This is a sim- 
plified notation, of course. The actual form of the integral 
operator can be found in (33l . l3~ij or within the literature 
on the Discrete Dipole Method [35| . Some theoretical 
results of importance to our discussion have been accu- 
mulated over the years. In this section we briefly review 
these results with the emphasis on understanding and 
physical meaning rather than mathematical rigour. 

B. Existence and uniqueness 

In operator notation equation ^ can be written as 

[I - GX]u = u in , (3) 

where / is the identity operator, G is the integral oper- 
ator with the Green's tensor kernel, X is a "diagonal" 
operator of pointwisc multiplication with the contrast 
function, u is the unknown total field, and u m is the 
known incident field. In our formulation the linear oper- 
ator I — GX belongs to the class of singular integral oper- 
ators [361 ] and its kernel is strongly singular as opposed to 
the weakly singular kernels of the corresponding integral 
operators in one- and two-dimensional scattering prob- 
lems. In [34[ the mathematical equivalence of the integral 
equation ((2]) to the frequency-domain Maxwell's equa- 
tions with the radiation boundary condition was shown 
(for Holder-continuous incident fields), and the necessary 
and sufficient condition was obtained for the existence of 
a solution. In the isotropic case with Holder-continuous 
contrast functions the solution of ([3]) exists if and only if 

e(x,w)^0, xel 3 . (4) 

Two sufficient conditions that guarantee the uniqueness 
of the solution are known. One is the presence of non-zero 
losses, i.e. Ime(x,w) > [34], EjjJ- I n the lossless case, 
Ime(x,w) = 0, a three times continuously differentiable 
permittivity function (with respect to all coordinates and 
in M 3 ) is also sufficient for the uniqueness [34| . 

Some authors prefer to work with an integro- 
differential form of equation ([2]), where the kernel of 
the integral operator is kept weakly singular (three- 
dimensional scalar Green's function of the Helmholtz 
equation), and the two spatial derivatives (grad-div op- 
erator) are kept outside the integral [H, di| • Although, 
the analysis is more complicated in that case and involves 
Sobolev rather than Hilbert spaces, in a condition 
similar to the existence condition (UJ) was also obtained. 



C. Spectrum 

One of the advantages of considering the operator 
/ — GX in its strongly singular form is that it nat- 
urally acts on the Hilbert space and one can apply a 
well established theory (3(| to analyze its spectrum. The 
physical importance of the eigenfunctions and eigenval- 
ues of / — GX stems from the fact that they describe 
the spatial spectrum of the field, similar to the eigen- 
modes of a closed resonator or the plane waves in the 
one-dimensional case. Detailed analysis of this problem 
was presented in [4lj | . where the spectrum was found to 
contain not only the usual eigenvalues, but a nontriv- 
ial essential part as well. This is also a purely three- 
dimensional phenomenon, related to the strong singu- 
larity of the kernel and not present in one- and two- 
dimensional scattering (42J. The difference between the 
eigenvalues and the essential spectrum can be explained 
as follows. Eigenvalues A and eigenfunctions u\ satisfy 

[/ - GX]u\ = Aw a, (5) 

where u\ belongs to the functional space in question - 
Hilbert space here. That is to say that the eigenfunc- 
tions may be viewed as some well-defined and, in prin- 
ciple, realizable spatial distributions of the field on D - 
an equivalent of the resonator modes. The exact location 
of eigenvalues in the complex plane is not known, only a 
bound is available: 

Ime(x, uj) — Ime(x, w)ReA+ 

(Rce(x, uj) - e ) ImA < 0, (6) 
|A| < oo. 

The last inequality follows from the boundedness of the 
operator. The distribution of eigenvalues inside this 
wedge-shaped bound depends on the scatterer and ap- 
plied frequency to. The essential spectrum, on the other 
hand, satisfies the Weyl definition: 

lim ||[/-GX]* n -A CBS *„|| =0, 

= 1, 

where, while each \I/ n belongs to the Hilbert space, the 
sequence {^ n } does not have a convergent subsequence, 
and does not converge to any function in the usual mean- 
ing of the word. Thus, the "essential" mode does not 
represent any physically realizable field distribution. It 
was shown in [4l[ that the essential spectrum of / — GX 
is given explicitly as 

A ess = £^, xel 3 , (8) 

i.e. it contains all values of the relative permittivity. 
While eigenvalues are "discrete", i.e. isolated points in 
the complex plane, the essential spectrum is obviously 
a dense set, if e(x, ui) is Holder-continuous as presumed 
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in the analysis. It was also shown in [43[ that the cor- 
responding sequence {^ n } is, in fact, a mollifier of the 
square root of the three-dimensional Dirac delta- function 
- a very exotic distribution localized around the position 
x c determined by the corresponding point of the essential 
spectrum A oss — e(x c ,o;)/£o- It is clear that definition 
(|7|) encompasses eigenvalues as well, since with eigenval- 
ues the sequence {^n} simply converges to a function it a 
from ([5]), which does belong to the Hilbert space. 

A connection of the essential spectrum with physics 
was suggested in [43| via the notion of the pseudospec- 
trum and its pseudomodes |44|. If the eigenvalues of 
/ — GX are defined as complex numbers A for which 



[i-gx -xir 1 



(9) 



then the e-pseudospectrum of I — GX may be defined as 
a set of all A ps satisfying: 



[I - GX - XpJ}- 1 



1 

> 

e 



(10) 



for some e > 0. By analogy, we extend the Weyl defini- 
tion as 

lim l|[7-GX]*„-A ps *„|| <e. (11) 

n— >N(e) 

In this case, if N(e) < oo, the sequence {^f n } will stop 
at some highly localized function *Bn, still in the Hilbert 
space and a physically realizable field distribution. It is 
a pseudomode though, since it ceases to be a function for 
e — » and N — > oo. 

Both the eigenvalues and the essential spectrum play 
important roles in resonant phenomena [45j |. If either an 
eigenvalue or the point of essential spectrum gets close 
to the zero of the complex plane, then a resonance is 
observed and most of the electromagnetic energy will 
be accumulated in the corresponding eigenfunction or a 
pseudomode, which therefore will determine the spatial 
distribution of the total field on D. Since all cigcnfunc- 
tions and pseudomodes are rapidly decaying, if contin- 
ued outside D, the resonances will generally lead to an 
increase of the field strength inside the scatterer and a de- 
creased field outside D - something one expects and ob- 
serves. The difference between the eigenvalue-based and 
the essential-spectrum-based resonances is in their phys- 
ical origins. The distribution of eigenvalues is strongly 
influenced by the geometry of the scatterer and operat- 
ing frequency. To get an eigenvalue close to zero one 
needs a proper combination of the size and frequency, 
similarly to a half-wavelength condition in a dipole an- 
tenna. Whereas, to get the essential spectrum close to 
zero we only need to have the permittivity of the ob- 
ject close to zero at some arbitrary point in D, inde- 
pendently of the object size and geometry. This is the 
main difference between a material-based (microscopic) 
and a geometry-based resonances. A natural microscopic 
resonance occurs if the permittivity exhibits anomalous 
dispersion. Since from the outside both resonances may 



look the same (e.g. decreased transmission or reflection), 
we can, indeed, view a metamaterial or a photonic crystal 
as a composite scatterer which mimics a natural micro- 
scopic resonance by a geometry-based one. 

Apart from describing the dominant field distributions 
during resonances, the use of eigenfunctions is rather lim- 
ited. This has to do with the fact that they are rarely 
given explicitly, but also with another unfortunate prop- 
erty of the operator I — GX, its non- normality. Re- 
call that an operator is non-normal if it does not com- 
mute with its own adjoint, i.e., [I — GX]*[I — GX] ^ 
[I — GX][I — GX]*. Non- normal operators are not uni- 
tary diagonalizable, hence, they cannot be diagonalized 
using their eigenfunctions. This is not specific to three 
dimensions and seems to be a fundamental property of 
the frequency-domain electromagnetic scattering, with 
its true physical significance yet to be uncovered. 



D. The scattered field 

When nondestructive measurements are performed, 
one is usually able to measure the field scattered outside 
D only. The scattered field is defined as E sc = E - E in . 
Once the total field E(x, u>) on D is obtained by solv- 
ing the integral equation @, the scattered field on some 
measurement domain S is obtained by simply evaluating 
the integral: 



E sc (x,w)= J G(x-x',o;)x(x',w)E(x',w)dx', 



(12) 



x G S. 



In operator notation we shall write 



RXu : 



(13) 



The integral operator R is not singular, not even weakly, 
since x / x' in the argument of the Green's tensor. 
Hence, we have an integral operator with a smooth, ab- 
solutely integrable kernel. It maps between the Hilbert 
space of functions with spatial support on D - object 
domain - and the Hilbert space of functions with spatial 
support on S - data domain. It is a compact operator 
and as such is very different from the previously consid- 
ered operator I—GX. The properties of R are important 
in inverse scattering and have been discussed in the cor- 
responding literature [13, HH • Compact operators have 
zero as the only accumulation point of their eigenvalues 
and often have zero as an eigenvalue too. For example, 
consider a single data-point so that it sc is just a complex 
number representing the measured complex amplitude 
of a single Cartesian component of the scattered field. 
Then, the integral formula (fT2l) - (fl~3| can be written as 
an inner product of two vector- valued functions: 



(Xu,r) = {w,r}, 



(14) 
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where r is a vector- valued function representing the com- 
plex conjugate of the k-th row of the Green's tensor 
Gfc m (x — x', to) with fixed xeS and x' taking values in 
D. We see that any wn 7^ orthogonal to r in the sense 
of the inner product would produce zero scattered field. 
Analogously, if we measure the scattered field at any fi- 
nite number of points, and the operator R represents the 
so-called semi-discrete mapping, then it is always possi- 
ble to find functions wn ^ 0, such that Rwn = 0. Such 
functions, often called non-radiating sources, belong to 
the null-space of operator R, denoted as Af(R). The con- 
sequence of all this is that the operator R does not have 
a bounded inverse, and solution of the so-called inverse- 
source problem Rw — u sc with respect to w may be not 
unique. 

In Q the measurements were reported inside a meta- 
material and the fields were compared (matched) not 
only on S, but on D as well. Obviously, however, no 
measurements are possible all over D without destroying 
the metal particles or perturbing the field distribution 
on their surfaces. Therefore, parts of D will remain inac- 
cessible, and the null-space of R will, in general, remain 
nontrivial. On the other hand, the case where instead of 
actual measurements one simulates numerically the fields 
on D, corresponds to RX = GX, and may constitute a 
complete data-set discussed in the following section. 



III. ANALYSIS OF THE INVERSE 
SCATTERING PROBLEM 

A. Effective inversion 

Let there be an inhomogeneous object with permittiv- 
ity e(x, u>), occupying the spatial domain D, illuminated 
by the incident field u ln , and producing the scattered field 
u sc over some data domain S. Naturally, we shall assume 
that for this real-world object there exists a unique solu- 
tion of the forward scattering problem ([3]) . Formally this 
solution may be written as 

u=[I -GX]~V n . (15) 

Hence, the scattered field is obtained as 

u sc = RX[I -CX^u™. (16) 

This is the main equation of the inverse scattering theory, 
where one wants to find the constitutive parameters of 
the object - the diagonal operator X - from the knowl- 
edge of the incident and scattered fields. We can see 
from (frnj) that the inverse scattering problem is not only 
nonlinear but also almost certainly ill-posed, due to the 
presence of the operator R. 

Our goal is to find another object, an effective scat- 
terer, occupying the same spatial domain D, but having 
a different permittivity function £ G f (x, uj) ^ e(x, u), such 
that the application of the incident field u ln will produce 



the same scattered field u sc as the original object. Hence, 
what we want is, in fact, 

RX el [I ~ GX cf ]-V n = RX[I - GXpV 1 . (17) 

The effective permittivity function can be simpler than 
the original, e.g. homogeneous (constant) over D. It can 
also be the result of spatial averaging or smoothing de- 
scribed mathematically as an application of some linear 
integral operator A on the original permittivity function, 
i.e. X { = AE — I = E c { — I, where E is the diagonal 
operator of relative permittivity. The only thing which 
really matters here is that E c { ^ E and X c f ^ X. Hence, 
we would like the inverse scattering problem (|16|) to have 
at least and at most two different solutions, X and X e {. 

This problem was considered in [33, |48[ in a completely 
different context - as a fast algorithm to determine the 
effective permittivity of a buried object with the goal 
to discriminate landmines from stones and other targets. 
To distinguish (|17p from the standard inverse scattering 
problem (|16p . the former is called the effective inversion 
problem. Obviously, it represents a generalization of the 
S-parameter retrieval method. 

It is possible to reformulate the standard effective 
medium theory (EMT) along the same lines. In EMT 
and in the derivation of macroscopic permittivity an av- 
eraged, smoothed total field inside the object is intro- 
duced. Let us denote this procedure as Bu, so that the 
averaging of the field inside the original highly inhomo- 
geneous scatterer, characterized by X , is obtained as 

Bu = B[I - GX]- 1 ^. (18) 

The main conjecture of the EMT and macroscopic elec- 
trodynamics is that the same field will be observed inside 
a suitably averaged effective object X e { = AE — I, i.e. 

[I - GX cf ]-V n = B[I - GX]-W n . (19) 

Since the field-averaging operator B is also a compact 
integral operator with a smooth kernel, the similarity 
with p7p is obvious. The two problems become equiv- 
alent, if we B-average the left-hand side of (fT!))) as well. 
Hence, the fundamental problems of effective inversion, 
S-parameter retrieval method, effective medium theory, 
and macroscopic electrodynamics are all mathematically 
equivalent up to the actual form of compact operators 
A, B, and RX, if both the true and the effective fields 
are averaged in the latter two approaches. From now on 
we shall concentrate on the S-parameter approach (effec- 
tive inversion) as the only one where the i?X-operator is 
explicit. 

B. Non-existence 

The first question we should ask ourselves is: does 
an effective scatterer exist at all? To this end we re- 
call a well-known property of inverse scattering prob- 
lems. While the compact operator R does not have a 
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bounded inverse, the inverse scattering problem may still 
have a unique, though unstable, solution [39(. The known 
uniqueness conditions in inverse scattering theory are suf- 
ficient, not necessary, and usually describe a particular 
set of incident fields and a set of measured scattered fields 
which together guarantee that there is only one solution 
to (|f 6|) . Let us call such a set a complete data-set. For 
example, in the isotropic case this set can be chosen as 
follows: the far- field pattern of the scattered electric field 
for all angles of observation, all angles of propagation of 
the incident time harmonic plane wave, three linearly 
independent polarizations and a single fixed frequency 
[37L [39l ] . Another possibly complete data-set may occur 
when the field due to the original scatterer is simulated 
inside D, i.e. RX — GX. Although, we are not aware 
of a rigorous proof of uniqueness for this rather artifi- 
cial data-set. The complete boundary measurements all 
around D may also be sufficient [49], as well as partial 
boundary measurements with some natural restrictions 
on the type of contrast functions [13, HH (proofs are avail- 
able for the static and diffusive regimes only). In any 
case, we may safely conclude that no effective scatterer 
exists on a complete data-set, as only the true scatterer 
will give the exact data-match. 

Thus, an effective scatterer, which matches the fields 
exactly, can only be found on some subsets of a complete 
data-set, i.e. for partial illuminations and/or observa- 
tions. Due to the sufficient nature of known uniqueness 
conditions, however, we cannot generally pinpoint a sub- 
set required for an effective scatterer to exist. On the 
other hand, we can already conclude that, if an effec- 
tive scatterer exists for a certain subset of a complete 
data-set, then either it does not exist or has a different 
effective permittivity on the complementary subset. In- 
deed, if we could find the same effective scatterer on two 
or more complementary subsets, then the inverse scat- 
tering problem would not be unique on the complete 
data-set - a contradiction. This is the reason behind the 
dependence of the effective permittivity on the illumina- 
tion/observation conditions. It shows that such depen- 
dence is a fundamental property of all effective models, 
not just the metamaterial slabs. 



C. Non-uniqueness 

Suppose that we were able to find an exact effective 
scatterer on some subset of a complete data-set. Is this 
effective scatterer unique? Since we are not sure about 
the required size of the subset, let us consider a very 
simple situation, where we have only one incident field 
and a single data-point. Then, from (fT4|) and (TiTl) we 
obtain 

(A cf [/-GA cf ]-V n ,r) = (X[/-GX]-VV} =« sc - 

(20) 

Obviously, any possible non-uniqueness of X e f is a prop- 
erty of the effective model itself, not of the particular 



data-point u sc . Although such non-uniqueness may be 
generally anticipated, it is rather difficult to prove. This 
is because, unlike the inverse source problem discussed 
at the end of the previous section, or inverse scatter- 
ing in the Born approximation, the present problem is 
non-linear. The fact that the operator R has a non- 
trivial null-space does not yet prove anything. Indeed, 
consider an effective model with the contrast function of 
the form x( x , w) = Xcf(w)p(x), where p(x) is the spa- 
tial profile function (in operator notation we shall write 
X c { = XcfP)- For a homogeneous model this profile may 
be defined as: p(x) = 1, x £ D; p(x) = 0, x ^ D; 
and p(x) is Holder-continuous across the boundary of D. 
For this model the problem reduces to finding just one 
complex number - the effective permittivity e c f or the 
effective contrast Xef- If we now apply the Born approx- 
imation, i.e., 



( X efP[I ~ XefGP]- 1 ^"^) » Xcf<iVV) 



(21) 



then the value of the effective contrast is uniquely deter- 
mined by a single data-point and a single incident field, 
and is given by 



Xef 



(Pu h \r) 



(22) 



However, if we do not make any approximations, then 
Xef is, in general, non-unique. To show this we start by 
comparing the following eigenvalue problems: 



GPu\ = Xu Xl 
[I - XeiGP]u x = X x u x . 



(23) 
(24) 



Using the same set of eigenfunctions we deduce the rela- 
tion 



\ x = 1 - Xef A. 



(25) 



Consider incident fields with different spectral content. 
First, we assume that our external source creates the 
field, which looks like one of the eigenfunctions: 



a x u x . 



(26) 



Now, using (j2"0)) and (f2"4"l) - ([2"5)) . we equate the data from 
two scatterers, characterized by Xef and x' c f — a Xcf, re- 
spectively, 

aXet - a x (Pu x , r) = -^L^a x (Pu x , r). (27) 
1 - ax c f A 1 - Xef A v 1 

This leads to the equation 

a 1 



1-aXefA 1— Xef A' 



(28) 



which has only one solution, a — 1. Thus, with a single- 
mode incident field this effective model is unique, i.e. 
Xef is uniquely determined by a single data-point u sc . 
However, a realistic incident field, e.g. a plane wave or 
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the field of a dipole source, will contain many different 
modes u\. Consider, for example, the incident field of 
the form 

u m = ai«i + a 2 u 2 , (29) 

where both Ui and u 2 are eigenfunctions. Then, we ob- 
tain the following problem for a: 

aci t ac 2 ci c 2 

1 - aX[ 1 - aX' 2 ~ 1 - X[ 1 - X' 2 ' (30) 

where ci j2 = ai,-z{Pui t 2,r) and A' x 2 = x c fAi ;2 . This 
problem reduces to a quadratic equation with two roots: 

"i = 1, 

_ ci + c 2 — ciA' 2 - c 2 X[ (31) 
a2 ~ ciA 2 +c 2 A'i- Cl (A 2 ) 2 -c 2 (A'i) 2 ' 

showing that the effective model is non-unique. The loca- 
tion of the second solution depends on the exact balance 
of the modes in the incident field, location of the receiver, 
and the eigenvalues, which in their turn depend on the 
applied frequency and the geometry of the object. In gen- 
eral, there will be as many solutions as there are modes 
present in the incident field. And we may expect them 
all to be different, if the modes correspond to distinct 
eigenvalues. 

D. Approximate effective permittivity 

Above we have considered data from a homogeneous 
object and showed that the homogeneous model itself 
cannot always be uniquely inverted. One can expect, 
however, that addition of just a few incident fields and/or 
measurement points will cure that problem. Hence, it 
would seem a reasonable strategy to add, say, new re- 
ceivers when we deal with the data from an inhomoge- 
neous object as well, thus obtaining a unique value of 
the effective permittivity. However, here we run into the 
problem of non-existence discussed above. Although, an 
effective permittivity with a single incident field and a 
single data-point always exists (this can be proven using 
the technique of the previous subsection), it is different 
for different observation points and incident fields (re- 
call our discussion on complementary data-sets). Hence, 
in general, no single effective permittivity will fit all the 
data, even if it were just two data-points and a single 
incident field. 

This brings us back to the original purpose of the effec- 
tive permittivity as a simplifying approximation. While 
the exact effective permittivity appears not to exist, we 
can still find an approximate effective permittivity which 
minimizes the discrepancy between the scattered fields. 
We have to realize, however, that this discrepancy will 
generally grow as we add new data to the problem. From 
being exactly zero for many different values of e e f to some 
finite value at some possibly unique minimum, which 



corresponds to the approximate effective permittivity. 
Whether that minimum is definitely unique we do not 
know. 

Further, since the location of exact effective permittiv- 
ities, obtained for each separate data-point, depends on 
the eigenvalues of the scattering operator, which in their 
turn depend on the geometry of the object, we expect 
the value(s) of an approximate effective permittivity to 
be geometry-dependent as well. 



E. Singularities 

A popular way of enforcing uniqueness on the effective 
permittivity is by choosing a "physical" value of e c f as 
opposed to other, "non-physical" ones. However, when 
constructing an effective model all we should care about 
is the match in the scattered fields and the well-posedness 
of the effective model. From this point of view, some 
currently disregarded values of e e £ are perfectly accept- 
able. For example, effective permittivities with the nega- 
tive imaginary part are often dismissed as representing a 
"non-physical" medium with gain. First of all, pumped 
media with population inversion are no less physical than 
ordinary lossy materials, being, in fact, less exotic then 
the negative refraction media. Thus, the effective model 
with gain, if it matches the scattered field data, should 
not be dismissed simply because one does not like it. 
There are certain values of the effective permittivity, 
though, which make the effective model ill-posed and 
should be disregarded. These are the points of perfect 
resonances where the spectrum of the effective scattering 
operator / — GX c f contains A = 0. In (4f| it was shown 
that a perfect eigenvalue-based resonance does occur in 
media with gain. However, this happens only for some 
specific combinations of the permittivity, frequency, and 
geometry. This is what one strives for in the design of 
laser cavities. Hence, only certain "discrete" values of e c f 
with Ime c f < will cause the trouble. 

A question of specific importance to metamaterials is 
whether real negative values of e e f are acceptable. In 
PH it was argued that in this case an essential-spectrum- 
based perfect resonance may occur rendering the effective 
model ill-posed. Indeed, the essential spectrum is A cs = 
e(x, ui)/eo, x <E M 3 , and in the Holder-continuous case one 
needs to connect somehow the real unit with a negative 
real value. Thus, if we use the shortest path, it will 
pass trough the zero of the complex plane. However, in 
principle, it is possible to create or imagine a scatterer 
whose permittivity avoids the e(x, uj) = point on its 
way from the real unit to the chosen negative real value. 
For instance, an object with a thin lossy outer coating 
seems to be acceptable from that point of view. There is 
another problem however, which was not anticipated in 
[45| . It turns out that there is also a potentially infinite 
number of discrete real values of e e f spread from zero 
to negative infinity, which all cause perfect eigenvalue- 
based resonances and make the effective model ill-posed. 
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Indeed, consider an effective homogeneous scatter with 
the contrast Xef( x >k>) = Pi i- e - Ecf/^o = 2. Then, 

|| [/ - GP A/]- 1 1| = [ J - T^a Gp ] _1 II = °°> 

' (32) 

if A is an eigenvalue of / — GP. Therefore, for a general 
homogeneous effective scatterer with permittivity e e f we 
have 

- GX^W = \\[I -(scf/so-^GP]- 1 ^™, 

(33) 

if 

E 1 — A 

where A's depend only on the geometry of D and the 
applied frequency. Thus, no matter what kind of illumi- 
nation and/or scattering data we use in (JT7J) , the effective 
homogeneous model will be ill-posed for the above values 
of effective permittivity. 

Let us figure out where these values are on the complex 
plane. To do so we need to know the location of A's for a 
homogeneous object with permittivity e f/eo = 2. They 
will generally occupy the lower half of the complex plane, 
i.e. ImA < 0. Then, from (f34|) we conclude that the 
troublesome values of e G f / £o will also be in the lower half 
of the complex plane. Thus we confirm the anticipated 
problems for some discrete values of e e f exhibiting "gain" , 
i.e. Ime e f < 0. There are, however, also purely real A's. 
For example, in [52| . it was shown that in the quasi- 
static case cu — » 0, i.e., for objects much smaller than the 
wavelength of the incident field, all discrete eigenvalues 
will be concentrated inside the convex hull of the essential 
spectrum. Since the essential spectrum is now a segment 
of real line stretching between one and two (provided that 
the Holder-continuous p(x) is chosen correspondingly), 
that is where all discrete eigenvalues will be distributed 
in the quasi-static regime. From (|3"4")l we conclude that in 
that case there will be infinitely many discrete real values 
of £of/so spread all over the interval (— oo,0], such that 
the effective model is ill-posed. We expect (see numerical 
examples in [ill . [4|| and here) that for our base object 
with e c f/eo = 2 real eigenvalues are present at higher 
frequencies as well, and that the above conclusion is not 
specific to small objects. 



F. Effective models of higher complexity 

In a recent review article [49] an explicit connection 
has been made between the inverse scattering (Calderon) 
problem of reconstructing an object from boundary mea- 
surements and the problem of invisibility, discussed in 
the context of metamaterials. Metamaterials are be- 
lieved to produce almost arbitrary anisotropic functions 



of effective permittivity and permeability. Since an- 
other well-known result in inverse scattering theory tells 
us that anisotropic objects cannot be uniquely recon- 
structed from boundary measurements [491 ] , the question 
arises about the possibility of employing effective models 
that in some aspects are more complex than the original 
scatterer. For example, a composite object consisting of 
isotropic parts could be modelled as an anisotropic one. 

Thus, in terms of the effective inversion problem (|17[) 
we are now interested in the possibility of X c f being a 
more complicated operator than the original X, so that 
a data-set which is complete for the original scatterer is 
incomplete for the effective model. If only the permit- 
tivity is considered to be anisotropic, then X e f becomes 
a three-diagonal multiplication operator. If a magnetic 
contrast is present as well, then the forward scattering 
problem will look like 
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(35) 



where Gki are known integral operators [34j, X c ^ m are 
three-diagonal contrast-operators, and e and h are the 
electric and magnetic fields, respectively. Uniqueness of 
the solution is guaranteed, if at least one of the constitu- 
tive parameters has losses or, in the lossless case, all com- 
ponents of all tensors are three times continuously differ- 
entiable functions of coordinates. The necessary and suf- 
ficient condition on the existence of the solution of ([35]) is 
also derived in [34| for Holder-continuous tensor-valued 
functions of permittivity and permeability, and requires 

3 3 

o k e km (x,uj)e m ^ 0, 

fc=l m=l 

3 3 (36 ) 

in 

k—1 m—1 

x e M 3 , 

where 9k, k = 1,2,3 are the Cartesian components of 
an arbitrary real vector of length one. Obviously, this 
condition will cause singularities in the inverse scatter- 
ing problem not just along lines and curves as in the 
isotropic case, but over whole areas of the complex plane. 
In particular, equation (|36[) seems to exclude the values 
of permittivity and permeability required for perfect in- 
visibility [49j| . 

When considering effective models of higher complex- 
ity we do not have to limit ourselves to anisotropy. Why 
not consider spatial dispersion or even extra spatial di- 
mensions? While these models may seem over the top, 
they simply illustrate the main problem with this ap- 
proach. We agree with the authors of [22], |23j, l3l| that, in 
general, effective models of higher complexity contradict 
the basic goal of effective modelling, which was to sim- 
plify the original scattering model. We also doubt that 
these models can give us any useful information about 
various exotic phenomena as suggested, for example, in 
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[53j. Indeed, if a higher complexity effective model ex- 
ists, then it must be, by definition, completely equiva- 
lent to the original composite, but mundane scatterer. 
Conversely, any electromagnetic effect which makes this 
exotic model identifiable, i.e., different from the original 
scatterer, will violate the existence condition on the ef- 
fective model. 

In fact, it is quite easy to come up with a higher 
complexity model, which seemingly shows physics, not 
present in the original object. For instance, we know 
that the value of effective permittivity for a homogeneous 
model generally varies with the incident field. Hence, we 
can say that e c f, which solves (fl7|) . depends on u ln . Al- 
though, multiplying the incident field by a constant does 
not change the value of e c f , adding another source to the 
existing one will cause a variation in e c f. Mathematically, 
this can be expressed as 



E c[ K + <]/£ d K], 



(37) 



showing that e c f is a nonlinear functional of the incident 
field. Thus a linearly reacting inhomogeneous object may 
be viewed, to a certain extent, as a nonlinear homoge- 
neous scatterer. 

Finally, the present isotropic single-frequency model, 
where the effective permittivity is allowed to vary with 
frequency, even if the permittivity of the original scat- 
terer does not, is, obviously, also a model of higher 
complexity, if considered over the range of frequencies. 
Whether addition of multiple frequencies to the partial 
illumination/observation data-set makes it a complete 
data-set is an open question. It probably does, since 
a considerable improvement in the reconstruction of spa- 
tially inhomogeneous Lorcntz-type dispersive media was 
observed in [54| with the addition of just a few frequen- 
cies. 



IV. FINDING EFFECTIVE PERMITTIVITY 

It is well-known that an analytical solution for an effec- 
tive 3D scatterer of general geometry D and/or general 
profile p(x) is not available. Numerical solution of the 
effective inversion problem (|f 7[) requires an inverse of a 
N x N matrix obtained after discretization of the inte- 
gral operator featured in equation J2]). This is practically 
impossible for an electrically large object, as N is in the 
order of millions. Yet, the numerical solution of the for- 
ward scattering problem §3§ for some particular u ln can 
be found with an iterative method. In the effective inver- 
sion problem, however, such solutions must be found for 
many values of the effective permittivity. Thus, if solu- 
tion takes M iterations and we need to find it for K values 
of e c f, then we would need to carry out M K iterations. 
Each iteration is computationally roughly equivalent to a 
matrix- vector product. Although the special symmetries 
of the integral operator G allow for a significant simpli- 



fication of this process via the FFT algorithm, we still 
have a serious computational bottleneck here. 

In principle, we could try to minimize the number of 
trial e e f's when searching for the one which minimizes 
the data-discrepancy. For example, we could apply a 
Newton-type minimization algorithm. However, this al- 
gorithm is good for finding the minimum of a single- 
minimum functional, and should not be used with multi- 
minima problems, as the one at hand here. Not to men- 
tion that we would certainly prefer to visualize the norm 
of the data-discrepancy for e c f in a certain range and 
confirm the analytical predictions made in this paper. 

An algorithm, which circumvents the computational 
bottleneck of the effective inversion problem, was pro- 
posed in (5f|. It exploits the invariance of the Krylov 
subspace constructed by the Arnoldi algorithm with re- 
spect to a shift by the parameter Xrf- Practically this 
means that the computationally expensive iterations will 
be carried out only once, say for Xef = lj and the re- 
sulting Arnoldi- vectors and the associated M x M upper 
Hessenberg matrix may be re-used with different values 
of x c f. Thus, instead of carrying out MK iterations with 
the problem of the order N x N, we have to perform M 
iterations with the problem of the order N x N and find 
solutions of K linear problems of the order M xM. Since 
typically M <C N, the total number of operations is now 
much smaller, and we have what is called a reduced- order 
algorithm. 

In this section we shall present numerical examples for 
the following three objects: a homogeneous object, an el- 
ementary cell of a dielectric photonic crystal, and a larger 
rectangular piece of this crystal. In all cases we shall visu- 
alize the discrepancy between the measured (here - sim- 
ulated) scattered field from a true object and the field 
from a homogeneous effective model of the same outer 
shape. Namely, we shall plot the log 10 of the following 
functional: 



1 ) RP 



(4-0 



GP 



u in || 2 



(38) 



|w sc || 2 



representing a two dimensional non-negative function of 
the real and imaginary parts of e/sq. In all our examples 
we focus on a single data-point case, i.e., u sc is just a 
complex number (both the amplitude and the phase of 
the field are presumed to be measured). Hence, the norm 
in (|38|) is simply the square of the absolute value of the 
discrepancy. If one wishes to know what would happen 
with more data-points included in the calculations, where 
the norm becomes a sum of the single data-point norms, 
then one simply has to add the corresponding functionals 



First we consider a homogeneous rectangular object of 
a resonant size depicted in Fig. 1 (a). The dimensions 
in terms of the free-space wavelength are indicated in 
the figure. The source of the incident field is an elec- 
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FIG. 1: Finding the effective permittivity of a homogeneous object: (a) - scattering configuration; (b) - surface of the data- 
discrepancy functional for Receiver 1; (c), (d) - data-discrepancy functionals as two-dimensional images for Receivers 1 and 
2. Non- uniqueness of the homogeneous model is clearly visible in (c). The presence of a "stationary" minimum (circle) means 
existence of the exact effective permittivity, which is simply the true permittivity of the homogeneous object e/eo — 4. 



trie point-dipole situated one wavelength away from one 
the object faces. The dipole is polarized as shown in the 
figure. We simulate a single Cartesian component (in- 
dicated by the orientation of the receivers in the figure) 
of the scattered field, one wavelength away from the two 
opposite sides of the object, simulating the transmission 
and the backscattering measurement setups. The rela- 
tive permittivity of the true object is e/eo = 4. The 
effective scatterer is also homogeneous and has the same 
outer shape. The discrepancy between the scattered field 
of the true object and the field scattered by homogeneous 
objects with other values of relative permittivity is shown 
as a surface in Fig. 1 (b), where the horizontal axes are 
the real and the imaginary parts of the trial relative per- 
mittivities, and the height of the surface gives the value 
of the discrepancy functional for Receiver 1 . To visualize 
this functional on a 100 x 100 grid we have computed the 
data-discrepancy for K = 10000 different values of the 
complex relative permittivity, which would not be possi- 
ble without a reduced-order algorithm discribed above. 

As expected from our theoretical analysis, the func- 



tional has many maxima and minima. Minima corre- 
spond to the match between the data, hence, giving the 
effective permittivities. Maxima, correspond to the sin- 
gularities of l[55]). which, as predicted by equation 
are situated in the lower half of the complex (e/eo)-planc 
and along the negative real axis. It is easier to analyze 
the shape of the functional on a two-dimensional image 
as the ones shown in Fig. 1 (c), (d), where the height of 
the functional is now represented by the brightness of the 
pixel. In these images bright spots correspond to singu- 
larities, and dark spots are the effective permittivities. 

Also shown is a contour around an area of the (e/erj)- 
plane where the forward scattering problem is solved with 
sufficient accuracy. The norm of the residual of the for- 
ward problem is smaller than 0.01 inside the contour 
(here, above the curve), i.e. we have beteer than 1% 
accuracy there. The variation of accuracy with permit- 
tivity is the result of a fixed number of Arnoldi iterations 
applied everywhere (here and below we use M = 60 itera- 
tions) . We achieve an almost machine precision for small 
contrasts, i.e. for e/erj in the neighborhood of one, but ar- 
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rive at progressively larger residuals for larger contrasts. 
Also, we cannot count on any good accuracy in the im- 
mediate neighborhood of discrete singularities and for all 
non-positive real e/eo, as the forward problem becomes 
numerically ill-conditioned there. Achieving acceptable 
accuracy for larger areas of the (e/eo)-plane comes at a 
cost of more Arnoldi iterations, and in our case is limited 
to M = 60 by the available computing resources. On the 
other hand, the depicted contour gives the norm of the 
residual of the forward problem, i.e. the total error in the 
solution all over the spatial doain D and for all Carte- 
sian components of the electric field. Whereas, what we 
should be concerned about is the error in the computed 
scattered field at the receiver location only. That error 
is probably much smaller than the norm of the resid- 
ual on D. Therefore, we may expect that the computed 
reduced-order functionals can be trusted way beyond the 
outlined contour. 

The different minima of the discrepancy functional, 
some in the "physical" part of the (e/eo)-plane, see 
Fig 1 (c), illustrate the non- uniqueness of the homo- 
geneous effective model. If we compare Fig. 1 (c) and 
Fig. 1 (d), which correspond to different locations of the 
receiver, then we notice that all singularities remain at 
the same points, since they depend only on the shape 
of the scatterer and the applied frequency. At the same 
time all but one minima have changed their locations. 
The one which did not move was, of course, the permit- 
tivity of the original homogeneous object e/so = 4. This 
gives us a way to retrieve the unique permittivity of a 
homogeneous scatterer from just two data-points via de- 
tection of a stationary minimum. 

Application of the homogeneous effective model to an 
inhomogeneous elementary cell, similar to those used in 
photonic crystals (56|, is depicted in Fig. 2. We have 
two cylindrical rods, both with the relative permittiv- 
ity e/eo — 4, and a A/3 gap between their centers, see 
Fig. 2 (a). We consider three different locations of the re- 
ceiver, indicated in Fig. 2-a as Receiver 1, Receiver 2, and 
Receiver 3, correspondingly. The distance of all receivers 
to the nearest faces of the scatterer is set to one wave- 
length. Both the true and the effective objects are smaller 
with respect to the wavelength than in the previous ex- 
ample (in the horizontal cross-section) . Hence, we see less 
singularities and minima in the images of Fig. 2 (b), (c), 
and (d). The other significant difference is the absence 
of a minimum, stationary with respect to the changes in 
the measurement setup. This illustrates the anticipated 
non-existence of the exact effective permittivity in the 
inhomogeneous case. 

Finally we consider a larger sample of a dielectric pho- 
tonic crystal, similar to the one investigated in [56j in 
relation to the phenomenon of negative refraction. In 
Fig. 3 (a) the scattering configuration is shown, where 
all three receivers are one wavelength away from the cor- 
responding faces of the rectangular effective object. This 
effective model has, in fact, the same dimensions as the 
one considered in Fig 1. The frequency is chosen in such 



a way that the cyliners in Fig. 3 (a) form a triangular 
photonic lattice in the horizontal cross-section with the 
lattice period equal to A/3. Thus, we are in the vicin- 
ity of the bandgap for this lattice, which might be the 
reason behind the large losses exhibited by one of the 
effective permittivity minima, see Fig. 3 (d). Otherwise 
there are still multiple minima for each separate receiver 
location, and no stationary minima. Also the minima 
of Fig. 2 and Fig. 3 are different, showing that the ef- 
fective permittivity depends on the geometry and size of 
the sample. Of course, it would be interesting to see, if 
there is a convergence in the location of the minima for 
progressively larger pieces of this photonic crystal. Un- 
fortunately, we were not able to investigate this question, 
due to the limitations imposed on us by the computing 
resources. 

The Arnoldi algorithm used in [55j and here is still a 
little too costly in terms of computer memory as it re- 
quires storage of M vectors of size A, which are used as a 
basis for the total field. The reason for this choice of a ro- 
bust but computationally expensive algorithm is the non- 
normal complex-symmetric system matrix of our prob- 
lem. A possible alternative based on the Pade via Lanc- 
zos process was proposed in [57| for the two-dimensional 
effective inversion problem. It requires storage of only 
three vectors of size N. A further generalization of the 
method, suitable for inhomogeneous scattering models 
and working with the Maxwell equations in their differ- 
ential form was proposed in 58]. 

We have presented here only a few examples which 
illustrate the main conclusions of the theoretical part 
of this paper. In addition, we have performed a large 
number of numerical experiments with different types of 
objects, looked at different field components, tried dif- 
ferent types of incident fields (plane waves and Gaussian 
beams), considered measurements in the far- field zone, 
etc. Every time we would get images quite similar to 
the ones presented here, with a large number of virtu- 
ally unpredictable minima. Except for the low frequen- 
cies, where within the considered range of permittivities 
we would typically get only one minimum, which was 
more or less stable with respect to the receiver location. 
When we considered a larger sample of a photonic crys- 
tal (7 rows, 12-13 cylindrical rods in each) at a very low 
frequency, then the observed minimum was in a close 
agreement with the effective permittivity predicted by 
the Maxwell-Garnet mixing formula. 



V. CONCLUSIONS 

We have presented a generalization of the S-parameter 
retrieval method to three-dimensional finite objects and 
arbitrary illumination and observation conditions. We 
view it as a special kind of inverse scattering problem 
- an effective inversion problem. Many, if not all, con- 
clusions of this paper equally apply to the S-parameter 
retrieval technique, the original low-frequency effective 
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FIG. 2: An attempt to find the effective permittivity of an inhomogeneous object: (a) - scattering configuration; (b), (c), 
and (d) - data-discrepancy functionals for Receiver 1, Receiver 2, and Receiver 3, correspondingiy. Notice the presence of two 
"physical" solutions in (b) and (d), demonstrating the non-uniqueness, and the absence of a "stationary" minimum in (b)-(d), 
demonstrating the non-existence of the exact effective permittivity for inhomogeneous objects. 



medium theory, and even to the derivation of the macro- 
scopic Maxwell equations, as we have shown here the 
mathematical equivalence of all these problems up to the 
form of averaging/scattering operators. Of course, anal- 
ysis in 3D is much more complicated and implicit than 
in ID, where an explicit analytical solution for a homo- 
geneous slab is available. Nevertheless, straightforward 
application of the spectral analysis and basic results from 
the inverse scattering theory showed that the general 3D 
case is similar to the well-studied ID slabs in many re- 
spects. 

The "exact" effective permittivity, which provides the 
exact match between the scattered field from an effective 
homogeneous object and the original inhomogeneous one, 
exists only on a limited set of incident fields and a lim- 
ited number of observation points (angles). In fact, we 
can only be sure about the existence of this exact ef- 



fective permittivity, either "physical" or not, when we 
have a single incident field and a single component of the 
scattered field observed at a single spatial location. We 
have shown that addition of just one more observation 
location may already cause the non-existence of the ex- 
act effective permittivity. On the other hand, addition of 
sources and/or receivers may lead to the uniqueness of an 
approximate effective permittivity (although we are not 
able to prove it yet). However, the more data we take 
into account the more approximate (less accurate) an ef- 
fective model becomes. Also, the approximate effective 
permittivity is geometry-dependent, making it difficult 
to view metamaterials and other strongly inhomogeneous 
composites as continuous media suitable for carving out 
arbitrarily-shaped optical devices. 

The exact effective permittivity, if it exists, is non- 
unique. As was shown here, this is due to the non- 
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FIG. 3: An attempt to find the effective permittivity of a larger photonic-crystal sample: (a) - scattering configuration; (b), 
(c), and (d) - data-discrepancy functionals for Receiver 1, 2, and 3, correspondingly. Although this crystal is built out of the 
elementary cells considered in Fig. 2, the effective permittivities are now different as they are influenced by the geometry of 
the sample, which has the same dimensions as the homogeneous object considered in Fig. 1. 



linearity of the effective inversion problem and the non- 
trivial null-space of the scattered-field operator. In the 
single data-point case, the number of additional exact ef- 
fective permittivities depends on the spatial spectrum of 
the incident field - the broader this spectrum, the more 
non-unique is the solution of the problem. The spatial 
spectrum here is the spectrum of the scattering oper- 
ator, rather than the usual plane waves. An incident 
plane wave has, in fact, a very broad spatial spectrum. 
The location of additional solutions in the complex plane 
is determined by the incident field and the spectrum of 
the scattering operator and is pretty arbitrary. Hence, it 
is not always possible to choose one of the solutions on 
some "physical" grounds. 

The effective inversion problem becomes singular for 
certain values of effective permittivity. Many of them lie 
in the "non-physical" area of the complex plane and show 
negative loss. As such they do not cause much trouble. 
However, real nonpositive values of the effective permit- 



tivity important for negative refraction and invisibility 
must be excluded as well, as for those values the solution 
of the forward scattering problem either does not exist, 
or is not unique, or both. Luckily, the location of sin- 
gularities depends only on the applied frequency and the 
outer shape of the object and does not depend on the 
illumination/observation conditions. 

Although we have expressed some scepticism about 
the usefulness of effective models that are more complex 
than the original scatterer, we consider this topic to be 
worth pursuing. For example, it is imperative to know 
if there exists a complete-data set, perhaps, larger than 
the boundary data of the Calderon problem, such that 
would guarantee the uniqueness of the solution of the 
inverse scattering problem for a general inhomogeneous 
anisotropic scatterer. Also, the formal difference between 
the effective inversion method and the effective medium 
(and macroscopic electrodynamics) approach, although 
small, needs to be considered in more detail. Even with 
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respect to the completely equivalent (double-averaged) 
EMT we do not know if there are conditions ensuring the 
uniqueness of the averaged inverse problem (|19|) . Thus, 
depending on the actual form of the averaging operator 
B, it may still turn out that the problems of the EMT 



and the macroscopic electrodynamics have exact solu- 
tions under the illumination conditions which preclude 
the existence of the solution of the effective inversion (S- 
parameter retrieval) problem. 
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